## Plot raw data
## These results are reported in Figure 2 and 3
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

data4=data3
mean=aggregate(data4$SO2_MASS..tons., by=list(data4$Year), FUN=mean)
names(mean)=c("Year","SO2(tons)")

mean$National=c(31.3042658730158, 29.6508928571428, 25.3117063492063,
                21.2175595238095, 15.7621031746031, 14.8224206349206) # Obtained from EPA in July 2020
mean$National=mean$National*25/35 #transfer the unit

setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "SO2_Emission_Trend.jpg")
jpeg(file=file.save.path,
     width=1000,height=550)
ggplot(mean, aes(x = Year))+
  geom_point( size = 2,  aes(y = `SO2(tons)`, shape="Sample SO2 Emissions (tons)") )+
  geom_line(aes(y = `SO2(tons)`, linetype="Sample SO2 Emissions (tons)"), size = 0.6) +
  geom_point( size = 2,  aes(y= National, shape="National SO2 Trend (ppb)")) +
  geom_line(aes(y= National,linetype="National SO2 Trend (ppb)"),size = 0.6) +
  scale_x_continuous(breaks=c(2013,2014,2015,2016,2017,2018),
                     labels=c("2013-14","2014-15","2015-16",
                              "2016-17", "2017-18","2018-19"))+
  scale_linetype_manual("",values=c("solid","dashed")) +
  scale_shape_manual("",values=c(1,2))+
  scale_y_continuous(name = "Sample Average Emissions (tons)",sec.axis = sec_axis(~.*35/25,name="National Trend (ppb)")
  )+
  annotate("text", x=2017,y=20, label="National Standard 75 ppb", size=8,color='black', alpha=0.8)+
  theme(legend.position="bottom",panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 26), 
        legend.key.size = unit(3,'line'),
        axis.text=element_text(size = 20), axis.title = element_text(size = 26))
dev.off()




mean=aggregate(data4$NOX_MASS..tons., by=list(data4$Year), FUN=mean)
names(mean)=c("Year","NOX(tons)")

mean$National=c(38.23168239, 37.95259434, 36.6125,
                34.8447327, 35.13207547, 35.08372642)  # Obtained from EPA in July 2020
mean$National=mean$National*16/39 #transfer the unit

setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "NOX_Emission_Trend.jpg")
jpeg(file=file.save.path,
     width=1000,height=550)
ggplot(mean, aes(x = Year))+
  geom_point( size = 2,  aes(y = `NOX(tons)`, shape="Sample NOX Emissions (tons)") )+
  geom_line(aes(y = `NOX(tons)`, linetype="Sample NOX Emissions (tons)"), size = 0.6) +
  geom_point( size = 2,  aes(y= National, shape="National NO2 Trend (ppb)")) +
  geom_line(aes(y= National,linetype="National NO2 Trend (ppb)"),size = 0.6) +
  scale_x_continuous(breaks=c(2013,2014,2015,2016,2017,2018),
                     labels=c("2013-14","2014-15","2015-16",
                              "2016-17", "2017-18","2018-19"))+
  scale_linetype_manual("",values=c("solid","dashed")) +
  scale_shape_manual("",values=c(1,2))+
  scale_y_continuous(name = "Sample Average Emissions (tons)",sec.axis = sec_axis(~.*39/16,name="National Trend (ppb)")
  )+
  annotate("text", x=2017,y=15.3, label="National Standard 100 ppb", size=8,color='black', alpha=0.8)+
  theme(legend.position="bottom",panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 26), 
        legend.key.size = unit(3,'line'),
        axis.text=element_text(size = 20), axis.title = element_text(size = 26))
dev.off()




data_AOD=subset(data4, !(is.na(AOD)))
mean=aggregate(data_AOD$AOD, by=list(data_AOD$Year), FUN=mean)
names(mean)=c("Year","AOD")

mean$National=c(8.877633967, 8.703930012, 8.419675481,
                7.616690795, 8.035327357, 8.163252221)
mean$National=mean$National*0.13/9 #transfer the unit

setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "PM_Emission_Trend.jpg")  # Obtained from EPA in July 2020
jpeg(file=file.save.path,
     width=1000,height=550)
ggplot(mean, aes(x = Year))+
  geom_point( size = 2,  aes(y = AOD, shape="Sample AOD") )+
  geom_line(aes(y = AOD, linetype="Sample AOD"), size = 0.6) +
  geom_point( size = 2,  aes(y= National, shape="National PM2.5 Trend (ug/m^3)")) +
  geom_line(aes(y= National,linetype="National PM2.5 Trend (ug/m^3)"),size = 0.6) +
  scale_x_continuous(breaks=c(2013,2014,2015,2016,2017,2018),
                     labels=c("2013-14","2014-15","2015-16",
                              "2016-17", "2017-18","2018-19"))+
  scale_linetype_manual("",values=c("solid","dashed")) +
  scale_shape_manual("",values=c(1,2))+
  scale_y_continuous(name = "Sample Average AOD",sec.axis = sec_axis(~.*9/0.13,breaks=c(7.0,7.5,8.0,8.5), label=c(7.0,7.5,8.0,8.5), name="National Trend (ug/m^3)")
  )+
  annotate("text", x=2017,y=0.123, label="National Standard 12 ug/m^3", size=8,color='black', alpha=0.8)+
  theme(legend.position="bottom",panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 26), 
        legend.key.size = unit(3,'line'),
        axis.text=element_text(size = 20), axis.title = element_text(size = 26))
dev.off()


## Power Plants Location
PP=unique(data4[c("Plant_Code","Latitude","Longitude")])
coordinates(PP)=c("Longitude","Latitude")
PP=st_as_sf(PP)

US=read_sf("Data/raw data/cb_2018_us_state_20m/cb_2018_us_state_20m.shp")
US=subset(US, !(NAME%in%c("Hawaii", "Puerto Rico", "Alaska" )))

st_crs(PP)= st_crs(US)

setwd("Code Outputs")
wd=getwd()
file.save.path <- file.path(wd, "MAP.jpg")
jpeg(file=file.save.path,
     width=3000,height=1800)
ggplot() +
  geom_sf(data = US, fill =NA, lwd=3)+
  geom_sf(data=PP, alpha=1, aes(shape="Coal-fired Power Plant"), lwd=5, stroke = 3,show.legend =F, color="blue", fill="red")+
  scale_shape_manual(values = c("Coal-fired Power Plant" = 21))+
  labs(shape="")+
  theme(legend.text = element_text(size = 30),
        legend.title = element_text(size = 40, face = "bold"),
        legend.key.size = unit(10,'line'), 
        panel.background=element_blank(),
        panel.border = element_blank(),
        axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank())
dev.off() 
